#include<iostream>
#include<eigen3/Eigen/Dense>
#include<eigen3/Eigen/Sparse>
#include<fstream>
#include<sstream>

#define size 119425

using namespace std;
using namespace Eigen;

typedef Triplet<double> T;
typedef SparseMatrix<double> SpMatrix;

void ReadSpMatrix(ifstream& in,SpMatrix& m);
void ReadRhs(ifstream& in,VectorXd& rhs);

int main()
{
    clock_t start,end;
    ifstream MatrixIn("spmat2.csv");
    ifstream RhsIn("rhs2.csv");
    ofstream out1("solv2LU.csv");
    SpMatrix m(size,size);
    VectorXd rhs(size),X(size),r(size);
    ReadSpMatrix(MatrixIn,m);
    ReadRhs(RhsIn,rhs);
    MatrixIn.close();
    RhsIn.close();
    cout << m.rows() <<endl;
    cout << m.cols() << endl;
    //cout << rhs;
    SparseLU<SpMatrix> solver;
    start = clock();
    solver.compute(m);
    X = solver.solve(rhs);
    end = clock();
    out1 << X;
    r = m*X-rhs;
    cout << "LU: "<< r.norm()/rhs.norm() << endl; 
    cout << "running time:" << (end-start)/CLOCKS_PER_SEC << endl;
}

void ReadSpMatrix(ifstream& in,SpMatrix& m)
{
    vector<T> triplelist;
    int i,j;
    double value;
    string line;
    string newline;
    string unit;
    while(getline(in,line))
    {
        istringstream sin(line);
        newline = "";
        while(getline(sin,unit,','))
        {
            newline += unit + " ";
        }
        istringstream Newsin(newline);
        Newsin >> i >> j >> value;
        triplelist.push_back(T(i-1,j-1,value));
    }
    m.setFromTriplets(triplelist.begin(),triplelist.end());
}

void ReadRhs(ifstream& in,VectorXd& rhs)
{
    for (int i=0;i<size;i++)
    {
        in >>  rhs(i);
    }
}